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The original motivation to build a quantum computer came from Feynman 1 1 1 who envisaged 
a machine capable of simulating generic quantum mechanical systems, a task that is believed to 
be intractable for classical computers. Such a machine would have a wide range of applications 
in the simulation of many-body quantum physics, including condensed matter physics, chemistry, 
and high energy physics. Part of Feynman 's challenge was met by Lloyd [2| who showed how to 
approximately decompose the time-evolution operator of interacting quantum particles into a short 
sequence of elementary gates, suitable for operation on a quantum computer. However, this left 
open the problem of how to simulate the equilibrium and static properties of quantum systems. This 
requires the preparation of ground and Gibbs states on a quantum computer. For classical systems, 
this problem is solved by the ubiquitous Metropolis algorithm [3], a method that basically acquired 
a monopoly for the simulation of interacting particles. Here, we demonstrate how to implement a 
quantum version of the Metropolis algorithm on a quantum computer. This algorithm permits to 
sample directly from the eigenstates of the Hamiltonian and thus evades the sign problem present in 
classical simulations. A small scale implementation of this algorithm can already be achieved with 
today's technology. 
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1 Introduction 



Since the early days of quantum mechanics, it has been clear there is a fundamental difficulty in studying 
many-body quantum systems: the configuration space - Hilbert space - of a collection of particles grows 
exponentially with the number of particles. Many of the important breakthroughs in quantum physics dur- 
ing the 20th century have resulted from efforts to address this problem, leading to fundamental theoretical 
and numerical methods to approximate solutions of the many-body Schrodinger equation. However, most 
of these methods are limited to weakly interacting particles; unfortunately, it is precisely when the interac- 
tions are strong that the most interesting physics arises. Notable examples include high-T c superconductors, 
electronic structure in large molecules, and quark confinement in quantum chromodynamics. 

The configuration-space explosion problem is not unique to quantum mechanics: the task of simulating 
interacting classical particles is challenging for the same reason. It was only with the advent of computers 
in the 1950's that a systematic way of simulating classical many-body systems was made possible. In 
their seminal paper [3| Metropolis et al. devised a general method to calculate the properties of any 
substance comprising individual molecules with classical statistics. This landmark paper is a cornerstone 
in the simulation of interacting systems and has had a huge influence on a wide variety of fields (see e.g. 
[4, 5 6 1). The Metropolis method can also be used to simulate certain quantum systems by a "quantum-to- 
classical map" [7]. Unfortunately, this quantum Monte Carlo method is only scalable when the mapping 
conserves the positivity of the statistical weights, and fails in the case of fermionic systems due to the 
infamous sign-problem. 

As the reality of quantum computers comes closer, it is crucial to revisit the original motivation of 
Feynman for building a quantum simulator and to develop a general method, suitable for quantum com- 
puting machines, to calculate the properties of any substance comprising interacting quantum molecules. 
Such an algorithm would have a multitude of applications. In quantum chemistry, it could be used to com- 
pute the electronic binding energy as a function of the coordinates of the nuclei, thus solving the central 
problem of interest. In condensed matter physics, it could e.g. be used to characterize the phase diagram 
of the Hubbard model as a function of filling factor, interaction strength, and temperature. Finally, it could 
conceivably be used to predict the mass of elementary particles, solving a central problem in high energy 
physics. 

The seminal work of Lloyd [2 1 demonstrated that a quantum computer can reproduce the dynamical 
evolution of any quantum many-body system. It did not address, however, the crucial problem of initial 
conditions: how to efficiently prepare the quantum computer in a state of physical interest such as a thermal 
or ground state. Ground states could in principle be prepared using the quantum phase estimation algorithm 
[8. 9 1 , but this method is in general not scalable, because it requires a variational state with a large overlap 
with the ground state. Methods are known for systems with frustration free interactions 1 10] or systems that 
are adiabatically connected to trivial Hamiltonians ifTTI . but such conditions are not generically satisfied. 
Terhal and Divincenzo [12| suggested two approaches of how a quantum computer could sample from 
the thermal state of a system. The first suggestion is also related to the metropolis rule, yet left open the 
problem of how one could get around the no-cloning result and could construct local updates which can be 
rejected. This shortcoming immediately leads to an exponential running time of the algorithm, as already 
discussed in their paper. The second approach of preparing thermal states is by simulating the system's 
interaction with a heat bath. However, this procedure seems to produce rather large errors when run on a 
quantum computer with finite resources, and a precise framework to describe these errors seems to be out of 
reach. Moreover, certain systems like polymers 1 13 1, binary mixtures [ 14| and critical spin chains |fT5l[16l 
experience extremely slow relaxation when put into interaction with a heat bath. The Metropolis dynamics 
solve this problem by allowing transformations that are not physically achievable, speeding up relaxation 
by many orders of magnitude and bridging the microscopic and relaxation time scales; this freedom is to a 
large extent responsible for the tremendous empirical success of the Metropolis method. 

In this paper we propose a direct quantum generalization of the classical Metropolis algorithm and show 
how one iteration of the algorithm can be implemented in polynomial time on a quantum computer. Our 
quantum algorithm is not affected by the aforementioned sign problem and can be used to prepare ground 
and thermal states of generic quantum many-body systems, bosonic and fermionic. Like the classical 
Metropolis algorithm, the quantum Metropolis algorithm is not expected to reach the ground state of an 
arbitrary Hamiltonian in polynomial time. The ability to prepare the ground state of a general Hamiltonian 
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in polynomial time would allow to solve QMA-complete problems. However, as a rule of thumb it always 
seems possible to define an update strategy for which the Metropolis algorithm thermalizes efficiently if 
the physical system thermalizes in polynomial time. There are no obvious reasons why the same should 
not be true for the quantum Metropolis algorithm. It also inherits all the flexibility and versatility of the 
classical method, leading, for instance, to a quantum generalization of simulated annealing J6). 

2 Summary of results 

In this section, we present a sketch of how the quantum Metropolis algorithm works. Details and general- 
izations will be worked out in later sections. 

To set the stage for the quantum Metropolis algorithm, let us first recall the classical version. We can 
assume for definiteness that the system is composed of n two-level particles, i.e., Ising spins. A lattice of 
100 spins has 2 100 different configurations, so it is inconceivable to average them all. The key insight of 
Metropolis et. al. was to set up a rapidly mixing Markov chain obeying detailed balance that samples from 
the configurations with the most significant probabilities. This can be achieved by randomly transforming 
an initial configuration to a new one (e.g. by flipping a randomly selected spin): if the energy of the new 
configuration is lower than the original, we retain the move, but if the energy is larger, we only retain 
the move with probability exp (f3(E id — E new )), where E is the energy of the configurations and (3 the 
inverse temperature. 

The challenge we address is to set up a similar process in the quantum case, i.e., to initiate an ergodic 
random walk on the eigenstates of a given quantum Hamiltonian with the appropriate Boltzmann weights. 
In analogy to a spin flip, the random walk can be realized by a random local unitary, and the move should 
be accepted or rejected following the Metropolis rule. There are, however, three obvious complications: 1) 
We do not know what the eigenvectors of the Hamiltonian are (this is precisely one of the problems that 
we want to solve). 2) Certain operations, such as energy measurements, are fundamentally irreversible in 
quantum mechanics, but the Metropolis method requires rejecting, hence undoing, certain transformations. 
3) One has to devise a criterion that proves that the fixed point of the quantum random walk is the Gibbs 
state. 

To address the first obstacle, we assume for simplicity that the Hamiltonian has non-degenerate com- 
mensurate eigenvalues Ei, and denote the corresponding eigenvectors \ipi). In the supplementary mate- 
rial, it is shown that those conditions are unnecessary. We can make use of the phase estimation algo- 
rithm ifTTl [T8l [8] [19] to prepare a random energy eigenstate and measure the energy of a given eigen- 
state. Then, each quantum Metropolis step (depicted in Fig. [TJ takes as input an energy eigenstate 
with known energy Ei, and applies a random local unitary transformation C, creating the superposition 
Cl^i) = Sfe^fcl 7 / 1 *:)- C cou ld be a bit-flip at a random location like in the classical setting, or some 
other simple transformation. The phase estimation algorithm is now used in a coherent way, producing 
Tlk x k\i } k) \Ek) ■ At this point, we could measure the second register to read out the energy Ek and accept 
or reject the move following the Metropolis prescription. However, such an energy measurement would 
involve an irreversible collapse of the wave function, which will make it impossible to return to the original 
configuration in the case of a reject step. 

Classically, we get around this second obstacle by keeping a copy of the original configuration in the 
computer's memory, so a rejected move can be easily undone. Unfortunately, this solution is ruled out in 
the quantum setting by the no-cloning theorem [20|. The key to the solution is to engineer a measurement 
that reveals as little information as possible about the new state, and therefore only slightly disturbs it. This 
can be achieved by a measurement that only reveals one bit of information — accept or reject the move — 
rather than a full energy measurement. The circuit that generates this binary measurement is shown at Fig. 
Q] It transforms the initial state \tpi) into 



where p k — min (1, exp (—f3(Ek — Ei))). The state can be seen as a coherent superposition of accept- 
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ing the update or rejecting it. The amplitudes xi \/~fl correspond exactly to the transition probabilities 
\ x k\ 2 fk °f tne Class i ca l Metropolis rule. The measurement is completed by measuring the last qubit in the 
computational basis. The outcome |1) will project the other registers in the state \ipf)- Upon obtaining 
this outcome, we can measure the second register to learn the new energy Ek and use the resulting energy 
eigenstate as input to the next Metropolis step. 

A measurement outcome |0) signals that the move must be rejected, so we must return to the input state 
As \ipt) is orthogonal to to we actually work in a simple 2-dimensional subspace, i.e. a qubit. In 
such a case, it is possible to go back to the initial state by an iterative scheme similar to the one employed by 
Marriott and Watrous in the context of quantum Merlin Arthur amplification ||2T1 . The circuit implementing 
this process is shown in Fig. [2] In essence, it repeatedly implements two binary measurements. The first 
is the one described in the previous paragraph. The second one, after a basis change, determines if the 
computer is in the eigenstate |?/>i) or not. A positive outcome to the latter measurement implies that we 
have returned to the input state, completing the rejection; in the case of a negative outcome, we repeat 
both measurements. Every sequence of these two measurements has a constant probability of achieving 
the rejection, so repeating recursively yields a success probability exponentially close to 1 . 

The quantum Metropolis algorithm can be used to generate a sequence of m states \4>j), j — 1, . . . , m 
that reproduce the statistical averages of the thermal state pa = e~^ H / Z for any observable X: 

-J2{<t> j \X\4> j ) = TrXp + 0(l/i/m). (1) 

3=1 

To show that the fixed point of the quantum random walk is the Gibbs state, we developed the theory 
of quantum detailed balance. Let be a complete basis of the physical Hilbert space and let {pi} be 

a probability distribution on this basis. Assume that a completely positive map £ obeys the condition 

y/PnPm (ll>i\£(\lpn){i>m\)\lpj) = y/PiPj |£(|V>j) (ipi I) I Ipn) ■ 

Then a = 'Yl li Pi\'4>i){4'i\ is a fixed point of £ . The quantum detailed balance condition only ensures that 
the thermal state pc is a possible fixed point of the quantum Metropolis algorithm. The uniqueness of 
this fixed point as well as the convergence rate to it depend on the choice of the set of random unitaries 
{C}. If the set of moves are chosen such that the map £ is ergodic, the uniqueness of the fixed point is 
ensured. This condition can be satisfied by choosing {C} to be a universal gate set 11221 . The Metropolis 
step obeys the quantum detailed balance condition, if the probability of applying a specific C is equal to 
the probability of applying its conjugate C\ This can be seen as the quantum analogue of the classical 
symmetry condition for the update probability. In some cases it even suffices to just apply the same local 
unitary C at every step of the algorithm (see Fig. 0J. In this case, the single unitary C has to be Hermitian 
and has to ensure ergodicity. The local unitary can be seen to induce 'non-local' transitions between the 
eigenstates because it is followed by a phase estimation procedure. 

Even though an implementation of this algorithm for full scale quantum many-body problems may 
be out of reach for todays technological means, we have presented an algorithm that is indeed scalable 
to system sizes that are interesting for actual physical simulations. A small scale implementation of the 
algorithm that can be achieved with present day technology is presented in the later sections. Moreover, 
a discussion is included that sketches the basic steps necessary for a simulation of some notoriously hard 
quantum many-body problems. Like in the classical setting the convergence rate and hence the runtime 
of the algorithm is dictated by the spectral gap of the stochastic map. The scaling of the gap depends on 
the respective problem Hamiltonian and the choice of updates {C}. Just as for the classical Metropolis 
algorithm, efficient thermalization is of course not expected for an arbitrary Hamiltonian. This would 
allow one to solve QMA-complete problems in polynomial time ll23l l24l l25l . It is however expected that 
the algorithm will thermalize if the physical system of interest thermalizes. The inverse gap of the quantum 
Metropolis map for the XX-chain in a transverse magnetic field at T = with a simple single spin flip 
update as shown in Fig. |4] This plot indicates that the gap scales like 0(1/N) with N the number of 
spins, even at criticality. To prove a polynomial scaling of the gap for more complex Hamiltonians remains 
a challenging open problem. Also, it is well known that the choice of updates {C} can have a dramatic 
impact on the convergence rate of the Markov chain in the classical setting. Finding good updates in the 
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quantum setting is a very interesting open question, although the above example suggests that the problem 
might be simpler in the quantum than in the classical case. The algorithm can be seen as a classical random 
walk on the eigenstates of the Hamiltonian. All samples are thus computed with respect to the actual 
eigenstates. This is why our method is suitable for the simulation of fermionic systems by exploiting the 
Jordan - Wigner transformation l29l as discussed in [ 30 1 . The fermionic sign problem is therefore not an 
issue for the quantum Metropolis algorithm. It is worth noting that an additional quadratic speedup might 
be achievable using the methods of ll26l l27l l28l . 

3 Description of the quantum Metropolis algorithm 

In this section, we provide a more elaborate description of the quantum Metropolis algorithm. The funda- 
mental building block is the quantum phase estimation algorithm (see section[5]i; throughout this section we 
assume that the phase estimation algorithm works perfectly, i.e. given an eigenstate \ipi) of the Hamiltonian 
H with energy Ei, we assume that the quantum phase estimation circuit $ implements the transformation 

|^)|0) -> \^)\Ei) 

where is encoded with r bits of precision. The fact that errors inevitably occur during quantum phase 
estimation will be dealt with in section|4] The algorithm runs through a number of steps 0..4 and, just as in 
the classical case, the total number of iterations of this procedure is related to the autocorrelation times of 
the underlying stochastic map. As analyzed in the next section, this procedure obeys the quantum detailed 
balance condition and hence allows to sample from the Gibbs state. The different steps are also depicted in 

Fig.GS 

Step 0: Initialize the quantum computer in a convenient state, e.g. 1 00 ... 0) . We need 4 quantum registers 
in total. The first one will encode the quantum states of the simulated system, while the other 3 registers 
are ancillas that will be traced out after every individual Metropolis step. The second register consists of 
r qubits and encodes the energy of the incoming quantum state with r bits of precision (bottom register 
in Fig. QJ). The third register is the one used to implement the quantum phase estimation algorithm, also 
with r qubits (top register QJ). The fourth register is a single qubit that will provide the randomness for 
accepting or rejecting the Metropolis step. 

Step 1: Re-initialize the three ancilla registers and implement the quantum phase estimation based circuit 
depicted in Fig. QJ followed by a measurement of the second register. This prepares an eigenstate \ipi) with 
energy Ej and associated energy register \Ei). The upper ancillas are left in the state |0) r as we assumed 
perfect phase estimation. The global state is now 

IV^)l^)|o)|o) 

Step 2: The next step is depicted in Fig. QJ). Assume that we have defined a set of unitaries C = {C} that 
can be implemented efficiently; those will correspond to the proposed moves or updates of the algorithm, 
just like one does for instance spin flips in the case of classical Monte Carlo. Just as in the classical case, the 
exact choice of this set of unitaries does not really matter as long as it is rich enough to generate all possible 
transitions; the convergence time will, however, depend on the particular choice of moves. The unitary C 
is drawn randomly from the set C according to some probability measure d/j,(C). It is only necessary that 
the probability of choosing a C is equal to the probability of choosing C', i.e. dfi(C) = dpk{C^), as this is 
dictated by the requirement that the process obeys detailed balance, cf. section |4~2l 
The new state can be written as a superposition of the eigenstates: 

k 

Implement the coherent quantum phase estimation step specified in Fig. QJ, which results in the state 
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J24m ^j2 x imm\E k )\o). 

k k 

Note that E k is only encoded with a precision of r bits, so that in practice there will be a lot of degeneracies. 
Finally, implement the unitary W(Ek,Ei) (Fig. [TJi) which is a one-qubit operation conditioned on the 
value of the 2 energy registers: 



W(E k ,Ei) = f ^ l ^jl ik -Jr^jg) (2) 

f ik = min(l,exp(-/3(£ fe -£;))). (3) 

The system is now in the state 



E 4 \Jfim\Ei)\E k ) |l> + £ 4^/1 - fftfa) \Ei)\E k )\0) . 

k k 

For later reference, the product of the three unitaries C, the phase estimation step, and W is called U (see 

Kg. at). 

Step 3: Measure the single ancilla qubit in the computational basis. A measurement outcome 1 corre- 
sponds to an acceptance of the move and collapses the state into 



E a! fcV/jiV'*>i^)i^*)ii) 



In the case of this accept move, we can next measure the third register which prepares a new eigenstate 
\ipk), and follow that by an inverse quantum phase estimation step. This leads to the state 



1^)1^)10) 1 1) 



with probability proportional to 



x kV fl 



. This state will be the input for the next step in the iteration of 
the Metropolis algorithm: go back to step 1 for this next iteration. Note that the sequence E — !• Qi — >• L 
depicted in Fig. [3]exactly corresponds to this sequence of gates. 

A measurement |0) in the single ancilla qubit signals a reject of the update. In this case, first apply the 
gate W , and then go to step 4. 

Step 4: Let us first define the Hermitian projectors Qo and Q\, made up of the gates defined in step 2 — 3 
including the measurement on the ancilla: 



Qo = I7 t (l<8l<8l8i|0)<0|)l7 
Qx = t/+(I(g)]I®I(g)|l)<l|)[/ 

Let us also define the Hermitian projectors Pq and Pi as 



Po = £ E \^ a )(^ a \<E)\E l )(E t \<E)I<E}l 

Pi = Yl E l^aXVal® \Ei){Ei\®I®I 

Here equality (or inequality) means that the first r bits of the energies do (not) coincide. This measurement 
P a can easily be implemented by a phase estimation step depicted in Fig. [IJ;. 
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The fourth step of the algorithm now consists of a sequence of measurements (see Fig. First we 
implement the von Neumann measurement defined by P a . If the outcome is Pi, then we managed to pre- 
pare a new eigenstate \ip a ) with the same energy as the initial one \ipi), and therefore succeeded in undoing 
the measurement. Go to step 1. If the outcome is Pq, we do the von Neumann measurement Q a . Inde- 
pendent of the outcome, we again measure P a , and if the outcome is Pi, we achieved our goal, otherwise 
we continue the recursion (see Fig. 01. It happens that the probability of failure decreases exponentially 
in the number of iterations (see section [3~Tl i , and therefore we have a very good probability of achieving 
our goal. In the rare occasion where we do not converge after a pre-specified number of steps, we abort the 
whole Monte Carlo simulation and start all over. 



This finishes the description of the steps in the algorithm. 



3.1 Running time of the rejection procedure: 

Let us discuss the convergence of the reject step more closely. As already explained, the algorithm should 
prepare a new state with the same energy as the original one Ei in the case of a reject move. As shown in 
Fig. |3] we will do this by repeating a sequence of two different binary measurements Pi and Qi. The recur- 
sion stops, whenever the measurement outcome Pi is obtained, where Pi is the projector on the subspace 
of energy P, ; . Note that it is crucial for the algorithm that the initially prepared state E\ipi) |0 2r+1 ) is an 
eigenstate of the projection Pi. This is indeed the case, even if we take into account the fluctuations in the 
quantum phase estimation step discussed in the next section: the error that is generated by the fluctuations 
of the pointer variable can be accounted for if we verify the equality of the energy in P only up to r < r 
bits of precision. This allows to enlarge the eigenspace of P with approximate energy Ei, encompassing 
the fluctuations of the pointer variable. 

Here we will calculate the expected running time. The probability of failure to reject the move, given 
that we start in some state in the energy Ei subspace, after n > 2 steps, is given by the probability of 
measuring Pq after n subsequent binary measurements. Note that the commutator [PqQ s Pq, PoQs'Pq] = 
for all s, s'. Therefore, see Fig. [3] the probability of failure can be cast into the form 



fail/ \ \ * 

pi W = 2^ 



Tr 



m=0 

XVil ® lO^KO 2 ^ 1 !) EQ P (PoQiP,)" 1 (PoQoPo) 1 



(PoQoPo)" m (P)QiP)) m PoQoE 



(4) 



The full expression can conveniently be summed up to a single term: 



p{ aU (n) = (M0 2r+i |PQoPo 



p (52QsPoQs) p> 



P o QoE\i>i)\0 



2r+l\ 



(5) 



We now make use of the Lemma d73l as stated in section|7]and choose a basis in which the projectors Pi 
and Qi are block diagonal. Note that we reuse the same two pointer registers at each phase estimation step 
in the algorithm. This means that even though a realistic phase estimation procedure does not necessarily 
act as a projective measurement on the physical subsystem, the binary measurements Pi and Qi are still 
projectors on the full circuit. Therefore Lemma (l73l l can still be employed, even for a realistic phase 
estimation procedure. Without loss of generality, we assume that the rank of rank(Pi) = p is smaller than 
the rank of Qi which is equal to half the dimension of the complete Hilbert space (note that Pi projects on 
a single energy subspace). Assume that the unitary Uj brings P and Q to this desired form. This allows us 
to rewrite © asp{ ail (n) = (^|(0 2r+1 |PC/jP fa u(n)U J E\ij l )\0 2r+1 ) with 



Dfaii(n) = 



( D {l~D){D 2 + (I - Df) n 
- V / D(I-D)(D 2 + (I — D) 2 Y 





-^D(l-D)(D 2 + (1-D) 2 y 
D 2 {D 2 + D) 2 ) n 





0^ 



1 
1 J 
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Here, D denotes a p-dimensional diagonal matrix with only positive entries. Note that the state 
UjE\ipi) |0 2r+1 ) has complete support on the projection operator Pi. That is, as we stated earlier, the state 
is an eigenstate of Pi . this means that it only acts on the first upper left block. If we denote by < d* < 1 
the diagonal entry of D that gives rise to the largest entry in the upper left block of the matrix -Df a ;i(n), we 
can bound 

p fail (n) < cf (1 - d*){d* 2 + (1 - d*f) n . (6) 

We observe, that the probability of failure decays exponentially in n, for a n-independent d* . Let us 
maximize this expression over all possible values of d*, in order to obtain an absolute upper bound to the 
failure probability. Defining x — d* 2 + (1 — d*) 2 = 1 — 2c?* (1 — d*), we see that this probability may be 
bounded by ^^x n . This expression is maximized by choosing x = for which we have 

1 / 1 \ n 1 
Pfaii(n) < — -y « — — — . (7) 
2(n+l) \1 + ±) 2e(n+l) 

Hence, choosing n = 0(1/ e) recursion steps is sufficient to reduce the probability of failure to below e. 
We have to choose this e in such a mannar, that the probability of failure during a complete cycle of the 
Metropolis algorithm is bounded by a small constant number. 



3.2 Running time of the quantum Metropolis algorithm 

Let us discuss the runtime scaling of the full Metropolis algorithm. In general, there are three types of error 
one has to deal with when we consider the the runtime scaling of the algorithm. 

First, we are dealing with a Markov chain and hence there is an associated mixing error e rmx . The 
mixing error of the Markov chain is defined with respect to trace norm distance, as \\£ mm ^ [p ] — a*\\i < 
e m i X . Here m m i X denotes the mixing time, i.e. the number of times the completely positive map has to 
be applied starting from an initial state p to be e m i X close to the steady state a* of the Markov chain. 
The mixing time is determined by the the gap A between the two largest eigenvalues in magnitude of the 
corresponding completely positive map. The trace norm is bounded by 0411 

\\e m {p]-o*}\i <C exp (l-A) m , (8) 

for a map that obeys quantum detailed balance, where C exp is some constant that typically scales exponen- 
tially in the system size. The runtime, or the mixing time, scales therefore as 

Just as for classical stochastic maps one needs to prove that the gap is bounded by a polynomial in the 
system size for each problem instance individually to ensure that the chain is rapidly mixing. It is generally 
believed, that to prove rapid mixing for a realistic Hamiltonian is hard. However, the convergence rate 
of the classical Metropolis algorithm is in practice favourable if the physical system thermalizes; this is 
because the Metropolis steps can mimic the actual physical thermalization procedure, albeit with the added 
flexibility of unphysical moves that make thermalization orders of magnitude faster. It is expected that the 
same will be true for the quantum Metropolis algorithm as well. 

The second type of imperfection relates to the fact, that the reject part of a local move cannot be imple- 
mented deterministically. However, we already showed, cf. 13.11 that this probability can be made arbitrary 
small by increasing the number of iterations in the reject move. For all realistic applications one would 
choose a fixed n* so that one only attempts to perform n < n* reject moves before discarding the sample. 
We want to achieve an overall success probability of preparing a valid sample that is bounded by some 
constant c. What do we mean by that? As already stated the Metropolis algorithm allows one to sample 
from the eigenstates \ipi) with a given probability pi ~ exp (—(3Ei). Since our reject procedure can only 
be implemented probabilistically we have to choose a fixed number of times n* we try to reject a pro- 
posed update. The probability of failure p&ii(n) of rejecting a proposed update after n steps is bounded by 
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Pfaii(?i) < 2e(w+i) 1 see ®' ^ or ^ algorithm to work, we want the algorithm to produce a sample after 
m m ix applications of the map £ with a probability that is larger than a constant c. Hence the probability of 
failure after m m ix steps should obey (1 — pm{n*)) mmix > c. This condition is met if we choose 

2e(l — c) 

This means, that we have to implement for each Metropolis step at most n* measurements Pi and Qi, 
before we discard the sample and start over again. Note that this is a very loose upper bound for the 
actual number of reject attempts, since the probability of failure actually decays actually exponentially in 
n, however, with some unknown constant that is ensured to be smaller than unity. 

The third error relates to the fact that we are implementing the algorithm on a quantum computer with finite 
resources, e.g. a finite register to store the energy eigenvalues in the phase estimation procedure. This leads 
to a modification of the completely positive map £, whose fixed point a* now deviates from the Gibbs state 
pc by 1 1 (7* — pc\\ i < £*■ This error will be discussed in section|4] 



4 Fixed point of the algorithm and influence of imperfections 

In the previous descriptions of the algorithm we only considered the idealized case when we are able to 
identify each eigenstate by its energy label. When this is the case, the algorithm can be interpreted as a 
classical Metropolis random walk where the configurations of the system are replaced by the eigenstates 
of a quantum Hamiltonian. However, this picture falls short if we consider the more realistic scenario of a 
Hamiltonian with degenerate energy subspaces. The rejection procedure ensures in this case only that we 
end up in the same energy subspace we started from. We therefore need to investigate the fixed point of the 
actual completely positive map that is generated by the circuit. We will see that the quantum Metropolis 
algorithm yields the exact Gibbs state as its fixed point, if the quantum phase estimation algorithm resolves 
the energies of all eigenstates exactly. This is obviously impossible for non integer eigenvalues as one 
would need infinitely many bits just to write down the energies in binary arithmetic. However, we will 
show that this is not a real problem. A polynomial resolution will yield samples that approximate the Gibbs 
state very well, if the Markov chain converges sufficiently fast. For the error analysis we will assume that 
the ergodicity condition is met, and that the problem Hamiltonian we are trying simulate is such that the 
Markov chain is rapidly mixing. To be precise, for the error analysis we assume that the Markov chain 
is trace-norm contracting, see section 14.31 We previously discussed the errors that arise due to the finite 
runtime of the algorithm in section [3~2"l and the error due to the indeterministic rejection scheme, cf. section 
13.11 In this section we consider the error that is related to the implementation of the algorithm. Due to the 
implementation on a quantum computer three types of error arise. 

1. Simulation errors. The quantum phase estimation algorithm requires implementing the dynamics 
U = e~ lHt generated by the system's Hamiltonian for various times t. This can only be done within 
a finite accuracy. 

2. Round-off errors. The quantum phase estimation algorithm represents the system's energy in binary 
arithmetic with r bits. This unavoidably implies that the energy is rounded off to r bits of accuracy. 

3. Phase estimation fluctuations. As seen in Eq. ( T7TT >, given an energy eigenstate of the system, the 
quantum phase estimation procedure outputs a random r-bit estimate of the corresponding energy. 
The output distribution is highly peaked around the true energy, but fluctuations are important and 
cannot be ignored. 

The first error is related to the fact that cxp(itH) has to be approximated by a Trotter-Suzuki unitary. 

This error can be ignored as long as the necessary effort in the simulation time Th to make this small, 
scales better than any power of l/ejy with e# being this simulation error [19|. This first source of error can 
be suppressed at polynomial cost. Another way to tackle this error is to adopt the analysis done in 11281 . 

The second type of error is not a problem on its own. Suppose that each eigenvalue of H is replaced by 
its closest ?'-bit approximation. The corresponding thermal state would differ from the exact one by factors 
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of exp(/32~ r ). By choosing r ^> log f3, this error can be made arbitrarily small. Note that the simulation 
cost grows exponentially with r, which implies that our Metropolis algorithm has complexity increasing 
linearly with /?. 

Interestingly, such a problem is already present in the classical Metropolis algorithm [ 32 ] , as one imple- 
ments the Markov chain on a computer with a floating point error. As a stochastic matrix is non-Hermitean, 
a tiny perturbation of the stochastic map (by introducing floating point arithmetic) could in principle change 
the eigenvectors drastically. However, nobody ever seems to have encountered such a problem; this might 
originate from the fact that the detailed balance condition ensures that the stochastic matrix is well behaved. 

The third type of error is more delicate and is intimately related to the second type. Indeed, it is not 
correct to suppose, as we did in the previous paragraph, that quantum phase estimation outputs the closest r- 
bit approximation to the energy of the eigenstate. Rather, it outputs a random energy distributed according 
to Eq. d7Tl ), sharply peaked around the exact energy. This distribution can be sharpened by employing a 
method developed in 11331 : the idea is to adjoin rj + 1 separate pointers, each comprising r qubits, and to 
perform quantum phase estimation T) times on the system using each of the first rj pointer systems in turn 
for the readout. Then the median of the results in the r\ pointers is computed in a coherent way and written 
into the (77 + l)th pointer. The probability that the median value deviates from the true energy by more than 
2~ r is less than 2~ v 11331 . Given an eigenstate of H, this leaves two possible phase estimation outcomes, 
corresponding to the r-bit energy values directly below and directly above the true energy. Hence, the high 
confidence phase estimation algorithm acts as 

hM0)-H^}(a<(IAI) + \\E i ])) + 0(e- r >) (11) 

where |a!j(|_-E,J)| 2 + |aj([~-E7j~|)| 2 = 1 and \_Ei\ and \Ei\ are the two closest r-bit approximations to 
Ei. Despite this improvement, it is not possible to make the outcome of the quantum phase estimation 
procedure deterministic. In the worst case where the exact energy for a given eigenstate falls exactly 
between two r-bit values, the two measurements outcomes will be equally likely. Thus, what we described 
in the main text as projectors onto energy bins are not truly von Neumann projective measurements, but 
rather correspond to generalized (positive operator valued measure, POVM) measurements on the system. 



Phase estimation unitary and POVM To understand this, let us start by writing out the full unitary $ of 
the standard quantum phase estimation procedure as defined in section The unitary acts on the TV-qubit 
register that stores the state of the simulated system and a single r-qubit ancilla register that is used to read 
out the phase information. We write 



2-1 2-1 



$=EE M ' 8 isXfi. where M * = E f( E i> x - v)\*i)Q>i\- (12) 

y—0 x—0 j — 1 



Note that the function 



is complex valued. The operators M^ =0 constitute the POVM generated on the system state by the phase 
estimation procedure. The label x of the POVM denotes the r-bit approximation to the energy generated 
by the phase estimation procedure, whereas y corresponds to the initial value of the ancilla register. The 
map $ is therefore the full unitary of the phase estimation procedure. Due to ( fT3l it becomes clear that the 
estimate x of the eigenvalue gets shifted by an amount of y, if the ancilla register is not initialized to 
y = Q. 



4.1 The completely positive map 

We now investigate the actual completely positive map (cp-map) generated by all unitaries and measure- 
ments in more detail. The full map can be understood as an initialization step denoted by E followed 
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by successive P and Q measurements, as discussed in section [3] and illustrated in Fig. [3] Note that the 
projectors Qi depend on the random unitary C. For each application of the map we draw a random unitary 
C from the set C — {C} according to the probability measure dp(C). We therefore have to average over 
the set C. The cp-map on the system is obtained by tracing out all ancilla registers. As shown in the pre- 
vious section |2 the error obtained by cutting the number of iterations in the reject case to n* can be made 
arbitrarily small; we can therefore approximate the full map as an infinite sum 

£[p] = j \r A [LQ l E(p®\tf T+l ){Q 2r+l \)EQ l lS] (14) 
+ Tr^ [PiQ E (p®\0 2r+1 )(0 2r+1 \) EQ^] 

OO 1 

+ E ^ A [P 1 Q Sft P ...P Q Sl PoQoE 

n—1 81 ...s n —0 

{p® |0 2l - +1 )(0 2r+1 |) EQ P Q S1 P . . . P Q Q an Pi] dp(C). 

The projective measurements P s and Q s are comprised of several individual operations. We adopt a new 
notation: an unmarked sum over the indices written as small Latin letters, e.g. ki,pi,... is taken to run 
over all 2 r integer values of the phase estimation ancilla register. The projectors can be written as 

Q s = E &m£m%C® ® |Pi><Pa| ® R s {kiM), (15) 

fel,fe2 Pl,P2 

P* = E E Mll ] M%l®\k{){kx\®\p x )&\®I, 
p i = E E M ^ M H ® ® |Pi><Pa| ® I- 

fcl=fe2 Pl>P2 

As before, we used the convention that the first register contains the physical state of the system. The 
second register of r-qubits corresponds to the register that stores the eigenvalue estimates of the first phase 
estimation, the third register is again used for phase estimation and the last register sets the single condition 
bit. The last matrix is defined as 

R s {k u k 2 ) =W(ki,k 2 ^\s)(s\W(ki,k 2 ), (16) 

with W defined in (0). Furthermore, the first operation in the circuit, that prepares an eigenstate and copies 
its energy eigenvalue to the lowest register, is denoted by 

= E E M l^ M H ® ®r k 2 ){h\ ® \ Pl )(p 2 \ ®I, (17) 

fcl ,&2 PI >P2 

where © r denotes an addition modulo 2 r . For notational purposes we introduced another operation 

i=E E M^+MgC^lfciXfcil^biXpsI®^,^). (18) 

fcl,fc2 Pl,P2 

A successful measurement of Q\ at the beginning of the circuit, Fig. 12 followed by the operation L corre- 
sponds to an acception of the Metropolis update and a further clean-up operation that becomes necessary, 
when considering a realistic phase estimation procedure. 

If we define new super-operators A[p] and B n ({s n })[p], the cp-map on the physical system can be 
written as 

oo 1 

E[ P ] =A[p}+B [p}+J2 E B n(Un})[p\- (19) 
n—1 si ...s n —0 



n 



Here A denotes the contribution to the cp-map that corresponds to the instance, where the suggested 
Metropolis move is accepted. Each of the B n correspond to a rejection of the update after n + 1 sub- 
sequent Q and P measurements. These superoperators can be expressed as follows: 



A[p] = E fMO) mm(l,e-^- fel )) M£C <M^ M° kl p M^M^M^ M&. 

(20) 



ki,k 2 <i,Pi,9i 

Furthermore, 



b o[p] = EE E [ (Ol&farJ&ikuhM (21) 

fcl h,ri d;pi,p 2 ;qi,q2 C 

Ml f MH ^ M h f M h CM l\ f M °i P M t f M £ M rl f M n f M t » 



and 



S„({s„})[p] = ^ ^ / d/i(C) g fcl ({s„},{Z n+ i},{r n+ i}) 



fcl d,{l n + 1 }-{r n+1 } 

D d ki ({l n+1 }) p D d ki \{r n+1 }) 



(22) 



The operators D and the scalar function g in the definition of B({s n }) n are given by 

g kl (KI, { W. K+i}) = (0|i? (fc!, n)i? Sl (fcx, ra) . . . (h, r n+1 ) (23) 

R s «(k 1 ,l n+1 )...R sl (k u l 2 )R (k 1 ,h)\0) 

and 

-EE M VKT^Mf^M^CMZr^MZ- 2 ^ . . . (24) 

W + l}^fcl {P2n} 

M^M^&Mf^M^CMl^Mt . 

"1 "1 tl tl /Cl Kl 

This concludes the description of the completely positive map corresponding to one iteration of the Metropo- 
lis algorithm. 



4.2 Fixed point of the ideal chain 

To be able to make statements about the fixed point of this quantum Markov chain, we introduce (see 
section [6]i a quantum generalization of the detailed balance concept. As for classical Markov chains, this 
criterion only ensures that the state with respect to which the chain is detailed balanced is a fixed point. 
However, it does not ensure that this fixed point is unique. The uniqueness follows from the ergodicity of 
the Markov chain [35 1 36 1 and thus depends in our case on the choice of updates {C}, which can be chosen 
depending on the problem Hamiltonian. A sufficient (but not necessary) condition for ergodicity can easily 
be obtained by enforcing {C} to form a universal gate set, as will be shown below. 

In section [6]it is shown that a quantum Markov chain obeys quantum detailed balance, if there exists a 
probability distribution {pi} and a complete set of orthonormal vectors {|V>i)} f° r which 

\/PnPm{fa \£[\fai){fan\]\lpj) = y/PiPj {fan | £ [| Ipj ) {fa \ ] | fax) • (25) 

This condition together with the ergodicity of the updates {C} ensures that the unique fixed point of the 
quantum Markov chain is 

2™ 

a = Y,Pi\fa){fa\- (26) 
»=i 
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We therefore would like to verify whether condition d25l l is satisfied when we choose the pi equal to 
the Boltzmann weights of H and the vectors equal to the eigenvectors 

The condition (f25l > is linear in the superoperators. We can therefore conclude that, when each of the 
summands A and all the B's in ([T9l individually satisfy this condition, the total cp-map £ is detailed 
balanced. 

The idealized case would be met if we could simulate a Hamiltonian H with eigenvalues Ei that are 
r-bit integer multiples of ^p, or if we had an infinitely large ancilla register for the phase estimation. In 
this case, the operators M E would reduce to simple projectors H E+p on the energy subspace labeled by 
E + p. Hence 

M^M% = S p , q U E+p . 

Note that the 5 p . q ensures that after each P and Q measurement the second ancilla register used for phase 
estimation is again completely disentangled and returns to its original value. 

Furthermore, in the special case when the eigenvalues of the Hamiltonian are non-degenerate the pro- 
jectors reduce to II ^ = \4>i)(ipi\- In this case it can be seen that the dynamics of the algorithm reduce 
to the standard classical Metropolis algorithm that is described by a classical stochastic matrix that can be 
computed as 

s y = W# [hWiOh/>i>. 

For this special case it is obvious that the detailed balance condition is met. 

Let us now turn to the more generic case, when the energy eigenvalues are degenerate. We investigate 
each of the contributions to the completely positive map ( fT9l ). 

The accept instance: We first investigate the accept instance described by the operator A[p\. 

A{p}=J2 f MC) min(l, e -^- £l >) Tl E2 CTl El P IL El tflL E2 . (27) 

The detailed balance criterion (f25j) forp; = ^e~^ Ei and reads 

^ e ^(E, +Ej )/2^ ilAUi){ ^ mm) = L e -P^+^)/^AA[\^ m ){rln\]\il>i). (28) 

Note that the chain of operators begins with a projector H El and ends with a projector H E2 . The detailed 
balance condition reads therefore 

1 -P{E i+Ej )/2 I Mc) min ^ e - 0(El - Ei) ^ S El ^J Ei , Ej (^\C\^)(^\C^\^ m ) (29) 
= I e -^+ £ ™)/ 2 J dy.(C) min^e-^-^)) S El , E J E ^ Ej (^\C\^ m )(^\C^\^}. 
Due to the fact that ^e~ f3El min (l, e ~^ Bi ~ E ^) = ^e^ l3E ' min (l, e ~l3(Ei-E t )y this re( j uce s to 

/ d(i(C) ^i\C\^){^\C^ m ) = I dfi(C) (^|C|^ m )(^|Ct|^), (30) 
Jc Jc 

where the energies of the eigenstates have to satisfy Ei = E m and Ei = Ej. 

One sees that ( |27] | is satisfied when the probability measure obeys 

dfi(C) = dfi(tf). (31) 

If we consider an implementation that only makes use of a single unitary C for every update, we have 
to ensure that this unitary is Hermitian, i.e. C = . This symmetry constraint on the measure can be 
seen as the quantum analogue of the fact, that we need to choose a symmetric update rule for the classical 
Metropolis scheme. 
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The reject instance: We now turn to the reject case described by the operators B n ({s n })[p] ■ The 
rejecting operators also simplify greatly when we consider the case of perfect phase estimation. After each 
phase estimation step the second register disentangles due to the S PltPl+1 , we get 

B n ({s n })[p] = J2 E 9E (KMM.K+i}) / dn(C) D% {{l n+1 }) P D%> ({r n+1 }) . (32) 
The chain of unitaries and measurement operators in the operator D $2M reduces to 

D° E ({WiH = u E c^u ln+1 cu E c^ . ..u^u h cu E , (33) 

where is the projector on to the orthogonal complement of energy subspace E. Note that the first and 
the last projector in each chain of operators is H E . Hence, all elements 

(V»l|Sn(K})[|^><^|]|^> 

vanish, if all energies are not equal E\ = Ei = Ej = E m . We can therefore disregard the probabilities pi 
on either side of the detailed balance equation d25l l. The detailed balance condition thus reads 

{MBn{{Sn})[mbP 3 \M m ) = (^|£„(K})[|VO<^IM). (34) 

It is important that the function g E {{s n }, {^n+i}> { r n+i}) (ED is real - Due to this fact and furthermore, 
since all the individual operators R S (E, k) are Hermitian, we may exchange the ordering of the indices 
{ln+i}, {r n +i}- That is, we may write 

9e ({s„}, {l n +i}, {r n +i}) = 9e ({s n }, {l n +i}, {r n +i})* (35) 
= (0| J R°(fc 1 ,Z 1 ) t J R Sl (fc 1 ,Z 2 ) t .i? s -(fc 1 ,Z Tl+1 ) + i2 s "(fci,^ +1 )+ . . . J R^(fc 1 ,r 2 )+i? (fc 1 ,r- 1 )t|0> 
= 9E{{s n },{r n+1 },{l n+1 }) 

Furthermore, since the individual projectors and are of course Hermitian, we may write 

mB n ({ Sn })[^)(^ m ) (36) 

= E 9e (K},{W,K+i}) / dp(c) 5 EltEl , Ej , Em {^\D° El ({i n+1 }) IViXVil^/ (K+i}) \i>m) 

= E 9E ({sn},K +1 },{l„ + i}) / dfi(C) 6 El ^ E] , Em (^\D Ei \{r n+1 }) \1> m )mD° Bl ({l„+i}) |Vi> 

= (^■|B n ({ Sn })[V' m )(^]|^)- 



The last equality in (|36l l is precisely due to the fact that we can reorder the indices as previously discussed 
and that we are dealing with projectors on the energy subspaces. 

As already said, a possible set of updates that will ensure ergodicity in general is given by choosing 
{C} equal to a universal gate set. So for instance the set of all possible single qubit unitaries augmented 
with the CNOT gate would suffice to ensure ergodicity for an arbitrary Hamiltonian. To show this, we 
make use of a result proved in [36|, Proposition 3. For completeness, we just repeat the part of the proof 
that is relevant to us. 



Primitive maps A completely positive map £ is called primitive if for all states p there exists a natural 
number m so that, 

£ m [p] > 0. (37) 

This means that £ m [p] has to be full rank for some m. All primitive maps are strongly irreducible,i.e. 
ergodic. That is, if £ is primitive the map has a unique eigenvalue \(£) with magnitude |A(£ )| = 1 and a 
unique fixed point a* > of full rank. 
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Proof: By contradiction: Assume that £ is primitive but not ergodic. This means that one of the follow- 
ing holds: (a) a* is not full rank; (b) There is another a* that corresponds to A = 1, i.e. the eigenvalue 
is degenerate; or (c) there exists another eigenvalue with |A'| = 1. If (a) holds the channel can not be 
primitive, since for all m we have £ rn [a*] = a* which is not full rank. Now, if (b) we will be able to define 
an e = [A maa; (((T*) _1 / 2 o'*(cr*)^ 1 / 2 )] _1 so that a* - ea* > is not full rank and we are back in case (a). 
Furthermore, if (a) and (b) do not hold but (c), the only other eigenvalues of magnitude 1 can only be a p 
-th root of unity for some finite natural number p. This implies, however, that assumtion (b) holds for the 
p-th power £ p , and thus (a) follows. 

With this Lemma at hand, it is straight forward to proof the uniqueness of the fixed point. All we need 
show is that the cp-map £ is primitive. 

Uniqueness of the Fixed point If we choose the set of all possible updates {C} equal to a set of universal 
gates, then the Metropolis Markov chain is ergodic for all finite (3 < oo. 

Proof: If £ denotes the map defined in (19[ , according to d37l i all we need to show is that there is an 
m such that for every \if>) and every p (ip\£ m [p]\i[;) > 0. Since p can always be written as a convex 
combination of rank 1 projectors it suffices to choose p — \ip)(ip\. Furthermore we observe that all B n 
defined in ( fT9] > are positive, i.e. 

^\B n {{ Sl })[p\\^)>Q, (38) 

since this expression can always be written as the trace over the product of positive semi-definite operators 
for any p and \ip), see <TT~4T >. We can therefore disregard the contributions from the B n and focus only on 
the accept instance A of the map £, since by virtue of d3~8l we have 

(il>\e m [\<p)(<p\]\il>) > (iMIvXvIM- 09) 

We can thus write 

(^\A m [\^\M = (40) 

P m 

dn(C 1 )...dn(C m ) Yl n™n(i,e-K Ei +i- E J)\&\n Em+1 c m ...c 1 n El \<p)\ 2 

J E t ...E m+1 i=l 

> e -P{E max -E min ) f d ^ Ci) d ^ Cm )F^{d, . . . C m ). 



Here E max and E m i n denote the largest and the smallest eigenvalues of the problem Hamiltonian H 
respectively, and we defined the integrant F as 

F^(c u ...c m )= Mn Sm+1 c m ...Cin Bl |^)| 2 . (4i) 

E±...E m + i 

Note that the prefactor e -P( E ™°-x-E min ) ^ Qes not van j sn f or a jj finite /?, Since the integrant F is non- 
negative, we only need to proove that F does not vanish. Since we are drawing the C\ . . . C m from a set 
of universal gates we can always find a finite m, by virtue of the Solovay - Kitaev theorem (37], so that 
there exists a sequence of gates Ci that ensures that there is a sufficiency large overlap between \ip) and 
C m . . . C\ \ip). That is for a given e m , there exists a sequence of m gates, so that 



|C m ...CiH| 2 



(ip\ILE m+1 C m ...CiU El \^) 



E 1 ...E„ 



> 1 - Cm, (42) 



where we inserted resolutions of the identity J2e Hgj. Hence, at least one of summands in d42b has to be 
non-zero and thus F^_ v is strictly positive and does not vanish. Therefore, there exists an integer m so that 
the integral in the last line of d40b is strictly positive. Since d40b acts as a lower bound to (ip\£ m [\ip) (<p\]\i/j) 
we can conclude that £ is primitive. 
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4.3 Error bounds and realistic phase estimation 



Let us next return to a more general Hamiltonian that has a realistic spectrum. As was discussed earlier, 
a realistic phase estimation procedure introduces errors not only due to the rounding of the energy values, 
but more importantly due to the fluctuations of the pointer variable. For a completely positive map with 
realistic phase estimation the detailed balance condition d25l l will not be met exactly, but we can show that 
the condition is satisfied approximately. This will be sufficient for our purposes. 

In order to bound this error we adopt a standard procedure also used for classical Markov chains |38|. 
Throughout this analysis we assume that the completely positive map is well behaved and is contract- 
ing. Whether this assumption is satisfied depends on the mixing properties of the problem we consider 
and on the choice of updates. Therefore, these properties have to be verified for every problem instance 
individually. A quantum Markov chain is trace - norm contracting if it satisfies 

||£[p-ff]||i<7Ji||p-<ii, (43) 

where the constant 771 < 1 is the smallest constant, so that this inequality holds |38|. The constant 771 is 
often referred to as the ergodicity coefficient. Note that the map is considered contracting only when the 
constant is strictly smaller than unity. It can occur, for some pathologically behaved maps, that this constant 
is not strictly smaller than unity even though the map is rapidly mixing. However, this can be cured by 
blocking several applications of the channel together, leading to a new constant smaller than unity [39|. 



Error bound The error e* between the exact fixed point a* of the map £ and the Gibbs state pa = 
4 exp (—(3H) can be bounded by 

Ik* - PG\\ < r^—. (44) 
I-771 

Here 771 < 1 is the ergodicity coefficient of £ and e S9 the error that arises due to a single application of the 

map on p G , i.e. \\£[p G ] -pcj||i < e ss ■ 



Proof: The error e* can be written as 

m 

\\<T*-p G \\ = lim \\£ m [ PG ]~ Pg\\i< lim Y,\\e k \p G ]-e k - 1 \pG]\\i (45) 

fe=l 

^ ,• fc— in or i ii \\£\Pg] - Pg\\i 

< lim > 77? \\£[PG ]-pG 1 = : • 

m— too * — ' I — 7h 

k=l n 

Thus we only need to bound the error that occurs when we apply the map £ to the Gibbs state p G once. 
In order to bound this error, we will make use of the fact that the completely positive map satisfies the 
detailed balance condition d25l > at least approximately. Let us discuss what it means to satisfy detailed 
balance approximately. 



Approximate detailed balance Suppose we are given a completely positive map £ and an orthonormal 
basis {IV'i}}. To each state we assign a Boltzmann weight of the form {pi = ^e~@ Ei }. If this cp-map 
does not precisely satisfy detailed balance, but only an approximate form such as 

S/P^m£\\*n)tym\]\ll>3) = VP^l(^rn\£U 3 )m\^n) (l+O^ 3 )) , (46) 

we can give the following bound on the error, measured in the trace - norm, that occurs upon a single 
application of the completely positive map. 

\\£{p G } - Pg\\i <C(e s s) (47) 
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Proof: Let us define p — J^. pi \ifji) (tpi \ . Then due to d46l l we have 



(fln\£\pa]\iM = ^2pi(ipi\£[\<Pi}(iPi\}\i>m) = (48) 

i 

Vp1p^(1 + 0(e S9 )) Tr [£[\^ m )(M] = Pm (1 + 0(e S9 )) S ml . 

So the application of £ yields £[pa] = Pa- Note that the state p~Q is still diagonal in the same basis as pc 
and both of the probabilities pi of p~a relate to the original probabilities via pi =pi(l + 0(e sg )). Since pa 
and p~a are both diagonal in the same basis, it is straightforward to compute that \\p~o — Pg\\i < 0(e S9 ). 



Let us now verify the approximate detailed balance condition d46| i of the completely positive map ( TT9b 
for a realistic spectrum of the Hamiltonian H . First let us consider the standard phase estimation procedure. 
Since the actual eigenvalues may have arbitrary real values, we may not assume that the individual Al% act 
as projectors on the system. Note that even the combination of M^M% is not Hermitian anymore when 
p q. This is precisely due to the fact that the function f(Ej ,k — p) < TT~3T > is complex valued. An additional 
phase is imprinted on the system state. At first sight this seems to hinder any form of detailed balance in the 
eigenbasis of the Hamiltonian. It turns out, however, that the total expression on either side of the detailed 
balance equation is still real. Note that is diagonal in the eigenbasis of H and assumes the form 

2™ 

M^Ml = J2 f(E, , k - p)*f(E 3 , k - q)\^)(^ | . (49) 
i=i 

Hence, the phases in f(Ej, k — p)* f(Ej, k — q) cancel up to a total phase factor f^^j^ , which is 
independent of both k and Ej. This allows us to write 

where now SP^ = S^ q . Let us have look at a segment of the chain of operators as they typically appear in 
the superoperators A or B ( fT9l . The typical sequences look like 

, , p iir(p 3 -pi) 

. . . Ml 3 M? 2 C M^Ml 1 ... -> ... 4^ r S p k 3P2 C S p k 2Pl . . . (51) 

K2 «2 fcl fcl g*5r-(P3— Pi) «2 Kl • ' 

This leads us to the conclusion that in each of the operator sequences the phases that arise due do to imper- 
fect phase procedure cancel. The first phase associated to po is due to the initialization, whereas the last 
phase associated with d is canceled due to the measurement. This gives an additional explanation of why it 
is necessary to reuse the same pointer register for the phase estimation procedure each time. However, this 
comes at a cost as the realistic phase estimation procedure doesn't naturally disentangle the pointer register 
used for the next phase estimation anymore. Hence, the initial state of the ancilla register for the next phase 
estimation step may be altered. So after subsequent measurements using the same register the distribution 
function of the pointer variable spreads. 

We now consider what happens in the case where we use the high confidence phase estimation based 
on the median - method [33 1. As already stated, this method allows us to perform phase estimation where 
the pointer variable fluctuates at most in the order of 2~ r . All other fluctuations are suppressed by a factor 
of 2~'' and will therefore be neglected in the following. According to ( fTTT i we can replace the function 
f(Ej, k — p) by its enhanced counterpart ctE (k — p), which acts as a binary amplitude for the two closest 
r-bit integers to the actual energy Ej . As discussed earlier, the phases that arise due to the imperfect phase 
estimation algorithm cancel, if for each of the -q phase estimations the corresponding registers are reused. 
We are therefore left again with operators acting on the physical system that are diagonal and have 
only real entries. We will thus regard the amplitudes oe^A; — p) as real from now on. We will therefore 



17 



write 

SI" = >>£,(*:- p)a Ej (k - qM)m (52) 



2' 



3 = 1 



Let us pause for a minute and have a closer look at the operators S^ q . As stated previously the Sjl q are 
diagonal in the Hamiltonians eigenbasis and have only real entries. Hence, these operators are Hermitian. 
Furthermore, since o? E , acts as a binary probability distribution on the two 5 = 2~ r closest integers to §-r, 
we see that for a fixed Ej and a fixed q, the only possible two values for k are 



E£ 
2tt 



q and k^ 



2vr 



Conversely, the operator S^ q has only support on the subspace spanned by the eigenvectors \ipj) whose 
energies lie in the interval 

Ej e [{k + q)- 2-' r ; (k + q) + 2~ r ] n [(k +p)- 2~ r ; (k + p) + 2~ r ] . 

This allows a further conclusion. For a fixed k and q the operator does not vanish only if 

p£ [q-2- r+x ;q + 2- r+1 ]. 

The interpretation is as follows: the operator S^ 9 implements the action of a phase estimation and its con- 
jugate on the system. If the ancilla register was initially in the state \q) the full phase estimation process 
does not disentangle the ancilla register afterwords, if we have performed in an intermediate operation. We 
have seen previously in the analysis for the idealized phase estimation procedure, see section l4~2l that the 
inverse phase estimation procedure returns the ancilla register to its original value \q). Since the pointer 
variable fluctuates now, this is not the case anymore and the pointer register remains entangled with the 
simulated system. However, since we perform an enhanced phase estimation procedure, the allowed values 
for the ancilla register are bounded by p^ = q ± 2~ r+1 . Thus even though S^ q is not a projector anymore, 
the previously discussed conditions suffice to ensure approximate detailed balance. 

Let us now verify the approximate detailed balance condition for each of the summands in (1191 . 

The accept instance: We analyze what happens in the accept case indicated by the operator A[p\. Due 
to the cancellation of the spurious phases (|5T1 i this operator has the form 

A[ P ] = £ ]T J dp{C) nun(l,e^^- fa )) SfrC S%° p S°f C^Sff. (53) 

fci,fc 2 d,pi,qi C 



We now want to verify whether the approximate detailed balance condition is met, when we choose again 
Pi = 4 l3E ' and \ipi) as the eigenstate of H. We choose a symmetric measure, i.e. dp(C^) = dp(C), and 
verify the approximate detailed balance condition d46| i. The left side of the equation reads 

le-^+^^^IAn^X^lll^) (54) 
= E E \e-" {Ei+E ^\dp{C) min(l,e-^-^) S^^^SlfCS^) 

ki,k 2 d,pi,qi C 

- E E §e"^+^)/ 2 f dp(C) min(l, e -W-^) (^\C\^)(^\C\^ m ) 



fci,fc2 d,pi,qi 

Q-E l (k 2 - d)a El (k 2 -pi)a Ei {k\ - Pi)a Ei (ki)a Em (k2 - d)a Em (k 2 - qi)a Ej (ki - q\)a E] (ki). 



We are free to relabel all the summation indices fci, fc2, d, . . . to match it with the other side of the equation. 
The sequence 

fc 2 = k\ + d ->■ ( Pl=q j +d , \ -+ h = k' + d -> d = 2 r - df (55) 
t Qi = Pi + « J 



18 



does exactly this. Note that since cue- (k + 2 r ) = ag. (fc) the constant 2 r in the last step can be dropped. If 
we now consider the worst case scenario of the fluctuations of (fci), we see that k\ deviates at most as 
much as ki « f£ ± 2~ r+1 . The same is also true for k-i and k' 2 ,k[ respectively. Hence we can conclude 




We can therefore establish, that 

^ e -P(E l+Ej) /2 {%/MAm{ ^ mm) = }_ e -fi(E l+ E m) ,2 { ^ mm){ ^ m (1 + 0(e)) (57) 

with e = P±f2- r which can be fully controlled by adjusting the relevant free parameters. 



The reject instance We now turn to the reject case. The operators change accordingly. We consider 
the detailed balance condition for each of the full B n ({s n })[p\. Note that due to the previously discussed 
phase cancellations the operators ({Z n +i}) as defined in d24b assume the form 

D d kl {Un+i}) = E E 5* a "C t Sf;;f 2 "- 1 C5»;- lW »- a C'^..5Sf 2 (7 t 5» Pl C5P l l0 . (58) 

The analysis of the reject case is very similar to the exact case. We make use of the fact that all the functions 
5fei {{s n }, {^n+i}, { r n+i}) and a,Ei — P) are rea l' anc l tnat we can re l a bel the indices like we did in the 
exact analysis. We have to establish that 

}_ e -^ + E 3 )/2 ( ^ {{sn}) m {ip , (59) 

= ie^^+^)/ 2 <^-|B„({ S „})[|V m ><^|]l^) (1 + O(e)) , 

up to some e, that will turn out to be e = n^j-f32~ r . We again start by considering the left side of (|59l and 
show that it will be equal to the right side up the specified e. 

~e-K* +E WMB n aa n })[\hMMl>rJ (60) 

= E E 9 kl ({Sn},{ln+l},{rn + l}) [ dn(C)^ (Bi+Ei)/2 

fel d;{/„ + 1 };{r„ +1 } C 

&j\D d J ({r„ +1 }) \i> m )mK ({Wl}) l^>- 

We will first exchange the index sets {r n+ i} and {l„+i}- This is possible since the function g kl is real and 
we follow the same analysis we already performed in the case of the idealized phase estimation. Now we 
turn to the sequence of the relabeling of the index set d, k±, h, r\, a\, b\, . . .. Note that a, and bi are part 
of the definition of Df i ({l n+ i}) and D% ({r n +i}) respectively d58l . The relabeling sequence that does 
what we want reads 

k l = k' 1 +d-> (61) 

P2n = q' 2n + d X ^ / l n +l = l'n+1 + d 



Vln-1 


= Q2n- 


x + d 


q_2n-\ 


= P'ln- 


.i + d 




= <zi + 




1 91 







I2n = P 2n + d J \ r n+1 = r n+1 + d 

a n+1 = b> +1 +d) ... ( l i= l \ + d I, J i'i '/i '» J, _. ,f = -r -,/'. 
b n+ i = a n+1 + dj ^ n = r x + a 

For these replacements to work, it is important to note that the operators R s (k% , lj) depend only on the dif- 
ferences, i.e. R s (k\—li). The sequence of replacements therefore leaves the function g kl ({s„}, {l n +i}, {tYh 
unchanged. However, since we do perform 2n phase estimation processes for each of the superoperators 
B„ ({s n })> the variable k\ in the last process may fluctuate in the order of n2~ r+1 , as was discussed earlier, 
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and we may no longer assume that the statistical weights on either side of the equation are equal. Hence 
we know that for the worst instance k% is 6 = ±n2~ r+1 close to either energy Ei , Ej , Ei , E m . We can 
therefore see, upon evaluating (1591 , that the detailed balance condition for each individual B n is met up to 

ane = n^/32- r . 

We observe that the e increases linearly in the number n of subsequent P and Q measurements we 
make to reject the proposed update. For all realistic applications, as discussed in section 13.21 one would 
choose a fixed n* so that one only would attempt to perform n < n* reject moves before discarding the 
sample. Since we want to achieve an overall success probability of preparing a valid sample that is lower 
bounded by a constant c, we have to choose n* > 2e (i-c) • Here m denotes the number of times we have 
to apply the map £ to be sufficiently close to the desired steady-state. This is related to the gap A of the 
map £, cf. section l3~2l Hence in the end we can give an error estimate for a single application of the map, 
which is of the order 



5 Implementation 

In this section we describe how to efficiently implement the quantum gates required by our algorithm on 
a quantum computer. As is now standard in the literature, we assume that we can implement single-qubit 
operations, measurements of the observables a a , and elementary two-qubit gates, such as the CNOT gate 
with unit cost. 

The first nontrivial operation required by our procedure is a means to simulate the unitary dynamics 
e -itH g enera ted by a /c-particle Hamiltonian H . We assume that H can be written as the sum of s terms, 
each of which is easy to simulate on a quantum computer. The best way to do this follows the method 
described by Berry et. al. \ 19 1 and by Childs BP : this procedure provides a simulation of the dynamics 
e~ using a quantum circuit of length Th, where 

T H = cs 2 t N (log, (iV)) 2 9 V lo s( s2 *o/^) , (63) 

and c is a constant, s denotes the number of summands in H, < t < to, en is the desired error, and 
log„ (N) is the function defined by 

log,(JV) =min{r| log^Af)}, 

where log^\-) is the rth iterated logarithm. Now, for a typical Hamiltonian encountered in condensed 
matter physics or quantum chemistry, the number of terms s scales as a polynomial with N, the number 
of particles. Thus the length Th of the circuit scales better than any power of 1/eh and is almost linear 
with to and scales slightly worse than a polynomial in N. Thus we can simulate e~ ltH for a length of time 
t ~ p(N) and to precision en ~ l/q(N) with an effort scaling polynomially with N, where p and q are 
polynomials. 

The next operation required by our algorithm is a method to measure the observable H. This can be 
done by making use of the quantum phase estimation lfl7l[T8ll . which is a discretization of von Neumann's 
prescription to measure a Hermitian observable. First adjoin an ancilla - the pointer - which is a continuous 
quantum variable initialized in the state |0), so that the system+pointer is initialized in the state 1^)10), 
where is the initial state of the system. Then evolve according to the new Hamiltonian K = H <g) p for 
a time t, so the evolution is given by 

2 N 

e -iwm = |^.)(^| ® e~ itE ^. (64) 

Supposing that is an eigenstate 1^-) of H we find that the system evolves to 

e-^|^)|0) = |^)|x = t£? i ). (65) 
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A measurement of the position of the pointer with sufficiently high accuracy will provide an approximation 
toE 3 . 

To carry out the above operation efficiently on a quantum computer we discretize the pointer using r 
qubits, replacing the continuous quantum variable with a 2 r -dimensional space, where the computational 
basis states | z) of the pointer represent the basis of momentum eigenstates of the original continuous quan- 
tum variable. The label z is the binary representation of the integers through 2 r — 1. In this representation 
the discretization of the momentum operator becomes 

p = E 2 —r~- m 

i=i 

With this normalization p\z) = -§r\z). Now the discretized Hamiltonian K = H g) p is a sum of terms 
involving at most k + 1 particles, if H is a fc-particle system. Thus we can simulate the dynamics of K 
using the method described above. 

In terms of the momentum eigenbasis the initial (discretized) state of the pointer is written 



2 r / 

This state can be prepared efficiently on quantum computer by first initializing the qubits of the pointer in 
the state |0) ■ • • |0) and applying an (inverse) quantum Fourier transform. The discretized evolution of the 
system+pointer now can be written 

1 T - 1 



e- itH ^)\x = Q) = —J2 e- iE "" 3r \il> j )\z). (68) 

z=0 



Performing an inverse quantum Fourier transform on the pointer leaves the system in the state \ipj) 
where 



Thus we find that 



where 



x=0 \ z=0 ) 
2 r -l 

!</>>= ^2 fiE^lx), (70) 

\m, X )\* = ±- -A \ (7i) 



4 r • 2 

^ sin 



which is strongly peaked near x = L"^rJ ■ To ensure that there are no overflow errors we need to choose 
t < u^jj-. (We assume here, for simplicity, that H > 0.) 

It is easy to see that actually performing the simulation of K for t = 1 using the method of lfl9ll requires 
a product of r simulations of the evolution according to ® f° r 1) 2, 2 2 , . . . , 2 r ~ 1 units of time, 
respectively. 

We write $ for the unitary operation representing the complete quantum phase estimation procedure. 
Using $ it is straightforward to describe a procedure to approximate a measurement of H: we adjoin r 
ancilla qubits and apply $ and then measure the ancilla qubits in the computational basis, approximately 
projecting the system into an eigenstate \ipj) of energy, and resulting in a string x which is an r bit approx- 
imation to the value Ej/\\H\\. 

Finally, let us briefly discuss how to implement the unitary gate W(Ek, Ei). This is a single qubit 
unitary conditioned on two energy registers. That this conditional unitary can be performed efficiently 
follows by observing that one can efficiently compute the angle 9 = arcsin(e 2-( 2 t e_ ^)) m to a scratchpad 
register, conditionally rotate the answer qubit by this angle, and uncompute 9. 



21 



6 Quantum detailed balance 



In this section we discuss the implications of Quantum detailed balance. We use detailed balance as a tool 
to ensure, that the constructed quantum Markov chain has the desired fixed point. 

Definition: Quantum detailed balance Let £ denote a completely positive map, and let a be a den- 
sity matrix, then the £ is said to obey detailed balance with respect to a if the induced map, £ a \p\ = 
£[a 1 / 2 p(j 1 / 2 ] is Hermitian with respect to the Hilbert-Schmidt scalar product. That is, the map has to 
satisfy Tr [p*^ [</>]] = Tr [£ CT [/j]t0] for all complex square matrices p and <f). 

If the completely positive map obeys detailed balance, we can immediately infer several properties. 
First of all, since £ can be obtained from £ a by a similarity transformation, £ must have a spectrum that is 
real. Furthermore a is guaranteed to be a fixed point of the completely positive map: 

Lemma: fixed point Let a be a state and £ [p] = ^2 A^pAj^ a completely positive map that satisfies the 
definition for Quantum detailed balance with respect to <r, then a is the steady state of £. 

Proof: Consider the two maps, £ a [p] = ^ A^a 1 ! 2 pa 1 ! 2 A\ and £*[p] = ^ a 1 ! 2 A^pA^a 1 ! 2 . By 
definition £ a [p] = £* [p] for all p. Then 

s[a] = £An = e; [i] = ^ 1/2 E A l A ^' 2 = °- 

We will now derive a simple criterion to verify whether a given channel is detailed balanced with respect 
to a specific state. Suppose the basis in which the density matrix is diagonal is known, then the detailed 
balance condition can be checked in a straightforward manner: 

Lemma: Detailed balance criterion Let {l^i)} be a complete basis of the physical Hilbert space and 
let {pi} be a probability distribution on this basis. Furthermore, assume that a completely positive map 
= E^pAt obeys 

\J PnPm 

then a = J2 i Pi\ipi) (ipi\ and £ obey the detailed balance condition. Therefore a is the fixed point of £. 

Proof: Let £ a be defined with respect to a — ^2iPi\iJi)(tpi\- We need to verify, whether £ a becomes 
Hermitian with respect to the Hilbert-Schmidt scalar product. One immediately sees that 

Tr[p f = E hi ^ nm VPnPm (tpj\£{\lp n )(tp m \)\tpi) 

ij;nm 

= E Pji$nm,/PW('lpm\£(\i>i}(i>j\)\i>n) 
ij;nm 

= Tr[£ a (p)Ul 



7 Binary measurements and pairs of subspaces 

The key technical reason why it is possible to implement the reject move in the quantum Metropolis algo- 
rithm is related to a very special normal form in which two (non-commuting) Hermitian projectors can be 
brought. 
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Lemma: Jordan 1875 Let Pi and Qi be two projectors of rank(Qi) = q and rank(Pi) = p on a Hilbert 
space H = C™ with p + q < n. We assume w.l.o.g., that q > p. Then there exists a basis of H in which Pi 
and Qi can be written in the form 



Pi 



Qi = 



D n 



On— p,p 
Op,n— p ™n—p,n—p 

( Dp y/Djj£~~ t ., 

y/D p (I p -D p ) I„ - D, n 



\ 







Here, D is a p x p diagonal matrix with real entries < d\ < 









®n-(q+p),n-(q+p) / 



<d p <l. 



(73) 



Proof: We can always choose a basis of H in which the projector Pi can be written as 



Pi 



I 



p 



n - 



>.n—p ^n—p.n—p 

In any basis, a general rank g projector Q\ can be written in the form 

Arin 



■Bn—p^q 



At 



Bl 



(74) 



(75) 



Here A pq and B n - p q are rectangular matrices over C. We require that Qi is a projector: Q\ = Qi leads 
to the constraint 



4t 4 
^pq^pq 



B 



Br*. — T 



n-p.q^^-p.q 



In 



(76) 



We can now choose to perform a singular value decomposition of A pq 
Ub^bVq. The projector can thus be written as 



U a ^aVI and B n 



U A 
U B 



H a vIVb^b 
Z b V b V a Za ^b^ b 



U 







Ul 



(77) 



Note that U a and [/g are p- and (n — p) -dimensional unitary matrices respectively. Therefore, the total 
block diagonal unitary U a ® Ub leaves the projector Pi invariant. If we turn to equation d76b . we see that 
upon inserting the singular value decomposition, the matrix V = V\ Vb must satisfy 



(78) 



Note that both S^S^ and I q — S^Es are diagonal matrices, which are according to (F78b similar. If we 
assume w.l.o.g., that the singular values are non-degenerate, we conclude that V can only be a permutation 
matrix. The degenerate case can be covered by a continuity argument. If we define D = E^E^ and apply 
the appropriate permutations to the remaining submatrices, we are left with the desired expression for Qi, 
To make the binary measurements complete, we have to choose the complementary projectors as Pq = 
I — Pi and Qq = I — Q\\ obviously, those complementary measurement projectors have a very similar 
structure to Pi and Q\, 



8 An experimental implementation 

It is possible to implement the quantum Metropolis algorithm with todays technology for a simple 2 qubit 
example system. Here, we will show how the different building blocks of the quantum Metropolis algorithm 
can be represented with simple quantum circuits. For this we need to consider a quantum computer of 5 
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qubits. Let us assume that we want to simulate the Gibbs states of the Heisenberg ferromagnet on 2 spin 
l/2s, i.e. 

H 2 = -- (of ®<T2 + of ®ol + <r'i ®cr|) , (79) 

which is certainly one of the most interesting Hamiltonians for 2 qubits. With the appropriate energy 
offset, this Hamiltonian has the spectrum {0, 2}, where the eigenvalue is threefold degenerate. This is 
very good news, as it means that an exact phase estimation algorithm can be set up with just a single (qu)bit 
of accuracy. Such a phase estimation requires simulating the Hamiltonian for a time t — ir/2. One sees 
that this unitary corresponds exactly to the SWAP gate. That is, 

U (|) = e~^ H2 = SWAP. (80) 

In the quantum Metropolis algorithm, we need to implement the controlled version of this SWAP, which 
is the Fredkin gate. In |j40l , it has been shown how this Fredkin gate can be implemented efficiently using 
optics. A related gate, the so-called Toffoli gate, was recently realized in the group of R. Blatt with an 
ion trap computer BP . The second gate to be implemented is the controlled Metropolis unitary W. The 
Metropolis unitary can be implemented with two controlled R y rotations: 

W(0p) = Ry{-O )cX Ry(6 ) c , (81) 

where we have made use of the standard single qubit unitary, R y (6) = exp(—i^-a v ) and wrote X — a x . 
The temperature can be controlled by the angle Bp. Comparison with the original Metropolis unitary (O 
shows that we have to set cos(6^) = e _/3 . The full circuit is depicted in Fig. [5] Note that this circuit 
can be simplified, if we regard the lowest qubit as a classical bit, which is determined by the first phase 
estimation. It is possible to condition the remainder of the circuit on the first phase estimation result. Then 
the controlled Metropolis unitary W can be implemented by a single CNOT operation. 

Let us briefly recall the necessary steps that are needed to implement the algorithm for this five-qubit 
example crcuit. Since the phase estimation procedure is exact, the algorithm simplifies greatly and all 
assumptions for the steps described in section[3]are met. We will recall the steps again in this paragraph, so 
that the section is sufficiently self-contained. We will, however, be less general and focus the description on 
the two-qubit Heisenberg Hamiltonian. The qubits that comprise the circuit are labeled according to Fig.|5J 
even though the order in which they are written corresponds to the notation used in the remaining part of 
the paper. This means that the first register, labeled by | ^23, contains the physical state of the system from 
which we sample. The second register |)i contains the value of the first phase estimation as indicated by 
the operation E (see Fig.|5t). Register number three is comprised of the fourth qubit |-) 4 and is used for 
the second phase estimation procedure which is part of the unitary U ( see Fig. [5j>). Finally, the fourth and 
last register is given by the accept/reject qubit number 5, i.e. |-)s. 

Step 0: Initialize the full circuit to the inital state 

\*h) = |oo) 23 |o)i|o) 4 |o) 5 . 

After the initialization continue with Step 1. 

Step 1: We currently are in a state that is of the form 

I^O) = |^)23|0)l|0) 4 |0) 5 , 

where 23 is some arbitrary two-qubit state stored in the second and third qubit. Apply the phase estima- 
tion map E as given in ( see Fig. and measure qubit number 1, 

4 

|^o) -> Etyo) =£<V<M |^)2 3 |^)i|0) 4 |0) 5 apply E 

i=l 

-> \-4>i) = \ijji)2z\E i )i\Q)i\Q)^. measure qubit 1 

Here the denote the eigenvectors of the Heisenberg Hamiltonian with the energies marked by Ei £ 
{0, 1}. Go to Step 2. 
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Step 2: We start by drawing a random unitary C with respect to a uniform probability distribution from 
the set of Pauli matrices {erf, erf, erf, erf} acting on either of the two qubits labeled by 2 and 3. We now 
apply the corresponding unitary U of Fig. to the state \ipi)- 

4 4 
|V> 2 ) = =^X fcl /fe^fe)23|^)l|^)4|l)5 + 5Za;HV / l-/^IV'fc)23|^)l|^)4|0) 5 . 

fc=l fe=l 

Here, the a^.; denote the matrix elements of C in the eigenbasis of and stems form the W matrix (fJJ, 
which is implemented by d8H and depicted in Fig.|5]3. 
Measure qubit number 5. 

accept: If the measurement outcome is 1 the corresponding state is proportional to 

4 

\fp+) oc ^x k i\/Jkt\4>k)23\Ei)i\E k ) 4 \l) 5 . 

k=l 

We then measure the second phase estimation register comprised of qubit number 4. With a probability 
proportional to \xki\ 2 fki the resulting state will collapse to an eigenstate of the Hamiltonian that is of the 
form 

|Vfc>23|-E!t)l|-Bfc>4|l)5. 

Then go to Step 4. 

reject: Otherwise, if the measurement outcome will be and the state is proportional to 

4 

IV'-) °c y^fctVl - 7ki\4>k)23\Ei)i\E k )i\0) 5 , 

k=l 

we have to start the rejection procedure. Go to Step 3. 

Step 3: We need to reject the proposed update. To this end we have to implement the measurement 
scheme as indicated in Fig. [3] The first thing we need to do is to apply the adjoined U* of the unitary in 
Fig-H}? to \ip—). We are left with the state 

|^ rej -) = C/t|^_). 

Starting from this state, we implement the following measurement scheme. 

Measure the projector P s (Fig. [5}:), 

where the outcome s = 1 corresponds to the case where the two energies agree and the outcome s = to 
the the case where two energies disagree. 

success s = 1 This outcome heralds that the energies coincide and that we successfully returned to the 
state prior to the proposed update C. Hence we have returned to a state in the energy subspace. Go to 
Step 4. 

failure s = We have failed to return to the original energy subspace. To unwind the state and to return 
to the original state we have to introduce a further binary projective measurement Q s . The measurement is 
related to the unitary U in the following manner. First we apply U again, then we measure qubit and 
finally we apply U>. Hence the Q s measurement reads 

Q s = I 23 ®Ii®Il4® |s)(s| 5 U. 

We now have to alternate the measurements Q s and P s . That is, we now repeatedly apply Q s , disregard 
the measurement outcome and apply P s , 

Qs -> Ps\ s=0 -^Qs^ Ps\ s=Q ■ ■ ■ 

until we measure the projector Pi once. The corresponding plan of action is given by Fig. [3] The result 
Pi indicates, that we have successfully returned to the original energy subspace. Go to Step 4. 
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Step 4: To finalize the single application of the Metropolis rule, we have to clean up the ancilla registers 
and prepare them for a subsequent application. The current state of the system is of the form 



|^fc>2 3 |£i>l|£fc>4|s>5. 



This state has to be mapped to 



|^fc>23|£i>l|£fc>4|s>5 -> |^fc>23|0>l|0) 4 |0) 5 . 



This prepares a new input state for the subsequent application of the Metropolis rule. We now have to 
return to Step 1 . 

This completes the description of the Metropolis algorithm. The number of times the sequence of the 
Steps 1-4 has to be repeated before a valid sample is prepared is related to the mixing time of the algorithm 
( see section [3T2l >. Note, that a different choice for the unitaries {C} in Step 2 is also possible. The only 
requirements are that the probability of applying C is equal to that of applying C\ and that the updates 
allow transitions between all the eigenstates. 

9 Simulation of quantum many-body systems 

It would go far beyond the scope of this paper to give a faithful account on only the most eminent appli- 
cations of the quantum Metropolis algorithm to the simulation of quantum many-body systems. We will 
therefore give only a brief sketch in this section on how we expect that the devised quantum algorithm will 
aid in the computation of static properties of some notoriously hard problems in quantum physics that have 
eluded direct computation for large system sizes by classical means. Such problems are for instance the 
determination of the phase diagram of the Hubbard model, the computation of binding energies of complex 
molecules in quantum chemistry, and the determination of the hadron masses in gauge theories. Common 
to these problems is that the particles are strongly interacting fermions and bosons. We expect that it is 
this class of problems where our algorithm will be able to give the strongest contributions. At this point 
we would like to point out, that the quantum Metropolis algorithm is not plagued by the notorious sign 
problem, because the algorithm allows one to sample directly in the eigenbasis of the Hamiltonian. This 
can be done irrespectively of whether the degrees of freedom are bosonic or fermionic. 
In order to implement the quantum Metropolis algorithm for a specific many-body Hamiltonian H, we 
need to be able to perform the phase estimation algorithm efficiently. The central subroutine that needs 
to be implemented is therefore the simulation of the time evolution for the Hamiltonian H <g> p, as was 
discussed previously in section [5] The simulation method described in 1 19 1 relies on the fact that we are 
able to decompose the Hamiltonian into a sum of local Hamiltonians hk with H = ^2 k hk that can be 
simulated by themselves on a quantum computer efficiently. A method to rephrase fermionic or bosonic 
degrees of freedom in terms of the quantum computational degrees of freedom , that is in terms of qubits, 
is therefore needed. Such a program was devised in [30 43, 44] and we merely give a brief overview here 
and refer the reader to the corresponding references. 

The Hubbard model: The Hubbard model ||42| is based on a tight binding approximation that describes 
electrons in a periodic potential confined to move only in the lowest Bloch band. The Hubbard Hamil- 
tonian consists of a hopping term and an interaction term written in form of fermionic creation c\ and 
annihilation operators that act on a lattice site i in a regular lattice of N sites. 



This Hamiltonian has to be expressed in terms of spin degrees of freedom in order to be implemented in 
the standard quantum circuit formulation. The interaction term can be seen to be implementable directly 
since the particle density rii jC operator acts only locally and is bosonic in nature. The implementation of 
the hopping term is a bit more challenging. Consider for simplicity the hopping term for a single electron 




(82) 



<i,j>,cr 
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spin only. This part can be expressed in terms of the Jordan-Wigner transformation, cf. Fig. |6j as 

* E \ (*f (®&+i«SK + . (83) 

<i,j> 

once a specific order of the N lattice sites has been chosen. As is shown in Fig. |6]the unitary evolution of 
each individual summand can be implemented with a cost that scales at most linearly with the total system 
size [43 44 1 . More general fermionic Hamiltonians can be implemented in a similar fashion. 

Quantum chemistry A central problem in Quantum chemistry is the determination of molecule proper- 
ties. The major challenge is the determination of the electron binding energies that need to be computed 
in dependence of the nuclei position. The general approach to this problem is to solve the approximate 
Hamiltonian of the electronic degrees of freedom that arises due to the Born-Oppenheimer approximation. 
In this approximation the nuclei positions are external parameters in the electronic Hamiltonian. The calcu- 
lation of the molecule properties relies on the fact that the electronic energy can be determined efficiently 
in dependence of the nuclei position. In their paper [45], Kassal and Aspuru-Guzik show how a quantum 
computer could be used to determine molecule properties at a time that is a constant multiple of the time 
needed to compute the molecular energy. The algorithm relies on a black box that computes the molecular 
energy for every configuration. The quantum Metropolis method can function as this black box algorithm, 
which was missing so far. For the Metropolis algorithm to work, one needs to implement the phase esti- 
mation procedure for the chemical Hamiltonian of interest. It is shown in 1461 . that the phase estimation 
procedure can be implemented efficiently for a general second quantized chemical Hamiltonian. 

Gauge theories The current most common non-perturbative approach to QCD is Wilson's lattice gauge 
theory flTl . which maps the problem to one of statistical mechanics, where the Euclidean action now 
assumes the role of a classical Hamilton function. It is therefore reasonable to assume, that lattice gauge 
theories would also be the method of choice for the quantum Metropolis algorithm. However, the algorithm 
relies on a Hamiltonian formulation of the problem. Such a formulation is given by Kogut and Susskind's 
[48 1 Hamiltonian formulation of lattice gauge theories in 3 + 1 dimensions. Here the 3-dimensional space 
is discretized and put on a cubic lattice, while time is left continuous. The fermions reside on the vertices 
of the lattice, while the gauge degrees of freedom are put on the links. The physical subspace is required to 
be annihilated by the generators of the gauge transformation, i.e. all physical states need to satisfy Gauss's 
law. 

It turns out however, that this approach seems to be very hard to implement on a quantum computer. This 
is due to the fact that each of the links carries a Hilbert space that is infinite dimensional, namely the space 
of all square integrable functions on the corresponding gauge group SU(N). A finite approximation to this 
Hilbert space therefore leads immediately to a breakdown of the underlying symmetry. 
A different formulation of gauge theories, that does not suffer from this problem, is therefore needed. Such 
a formulation is given in terms of quantum link models introduced by Horn |49|. Brower et al. showed 
that QCD and in general any SU (N) gauge theory can be expressed as a quantum link model ||50l . In the 
quantum link formulation the classical statistical mechanics problem is replaced by a problem formulated 
in terms of quantum statistical mechanics in which the classical Euclidean action is replaced by a quantum 
Hamiltonian. The central feature is that the corresponding Hilbert space of the gauge degrees of freedom at 
each link is now finite. It suffices that each link of a SU(N) link model carries a single finite representation 
of SU(2N). This is achieved by formulating the problem in 4 + 1 dimensions, where the four physical 
dimensions correspond to the actual physical Euclidean space time, while the fifth Euclidean dimension 
plays the role of an additional unphysical dimension. The 4-dimensional Euclidean space time is discretized 
and lives on a cubic lattice. Furthermore, it was shown by Brower et al. [50], that the continuum limit is 
obtained by sending the fifth unphysical Euclidean dimension to infinity, which corresponds to preparing 
the ground state of the lattice Hamiltonian. It can be seen, that the 4 + 1 dimensional link models are related 
to standard gauge theories in 4 dimensions via dimensional reduction ifSTI . 

The full Hilbert space of the SU (3) gauge theory can be written as the tensor product of a 20-dimensional 
Hilbert space for each link of the lattice and the finite dimensional fermionic Hilbert space that resembles 
the quarks. In contrast to the standard lattice gauge theories the configuration space of the quantum link 
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model resembles that of quantum spin models. 

The physical spectrum, and by that the Hadron masses, of the 4-dimensional theory can be obtained from 
computing the correlation functions in the Euclidean direction on the ground state of the 4-dimensional 
lattice Hamiltonian. 
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Figure 1 : Fig. (a) The first step of the quantum circuit: the input is an arbitrary state | ip) and two r-qubit reg- 
isters initialized to 1 0) r . Quantum phase estimation $ is applied to the state and the second register. The en- 
ergy value in this register is then copied to the first register by a sequence of cnot gates. An inverse quantum 
phase estimation is applied to the state and the second register . Fig. (b) The elementary step in the quantum 
circuit: the input is the eigenstate |?/>i) with energy register \Ei) and two registers initialized to |0) r and |0). 
The unitary C is then applied, followed by a quantum phase estimation step and the coherent Metropolis 
gateW. The state evolves as follows: |^)|^)|0)|0) ->■ C|^)|ff, )|0)|0) = E fc 4l^fe>l^>|0)|0) -> 
Ek4\^m)\E k )\0) -> J2k 4Vfl\^h)\Ei)\E k )\l) + Ek4V^7 k \^)\Ei)\E k )\0) with/* = 
min (1, exp (—f3(Ei — E k ))). Fig. (c) The binary measurement checks whether the energy of the state 
\tp) is the same as the energy of the original one \4>i}- This is done by using an extra register containing 
phase estimation ancillas, a step that checks whether the energy is equal to Ei or not, and finally an undoing 
of the phase estimation step that preserves coherence. 
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Figure 2: The circuit corresponds to a single application of the map £. The first step E prepares an 
eigenstate of the Hamiltonian, The second step Qi , measures whether we want to accept or reject the 
proposed update. In the "reject" case the complete quantum circuit comprises a sequence of measurements 
of the Hermitian projectors Qi and Pj. The recursion is aborted whenever the outcome Pi is obtained, 
which indicates that we have returned to a state with the same energy as the input. Because each itera- 
tion has a constant success probability, the overall probability of obtaining the outcome P x approaches 1 
exponentially with the number of iterations. 



30 




Figure 3: Given an input state \tp), we first perform phase estimation to collapse to an eigenstate with known 
energy E. This graph represents the plan of action conditioned on the different measurement outcomes of 
the binary P and Q measurements. Each node in the graph corresponds to an intermediate state in the 
algorithm. One iteration of the map is completed when we reach one of the final leafs labelled by either 
accept or reject. The sequence E — > Qi — > L corresponds to accepting the update, all other leafs to a 
rejection. 
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Figure 4: Inverse gap of the quantum Metropolis map at T = as a function of the number of spins in 
a chain with Hamiltonian T~L = J2k XkXu+i + YkYk+i + gZk- The update rule is a single-spin flip Xi; 
remarkably, this single gate is enough to ensure ergodicity. The observed linear scaling indicates that, at 
least in the case of ID spin chains with nearest - neighbor Hamiltonians, the quantum Metropolis algorithm 
converges in polynomial time. 
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Figure 5: Fig. (a) describes the first phase estimation step of the circuit. Since the phase estimation 
of the two-qubit Heisenberg Hamiltonian can be implemented exactly by the Fredkin gate, a single phase 
estimation operation is sufficient. In Fig . (b) the elementary unitary of the circuit is depicted. The angle 
of the controlled-controlled R y {0fj) rotation needs to be chosen such that cos{Qp) = e _/3 . The final 
measurement P is depicted in Fig. (c). The first phase estimation has to be followed by a measurement 
which verifies that the two phase estimation bits are equal. The phase estimation is then undone so that P 
is a Hermitian projector. 



31 



1 

2 



■exp(ieiTj) ■ 



-V 12 



exp(-ieof o-fof; 



-Kit 



13" 



ui- 



Figure 6: A fermionic many particle Hamiltonian can be simulated on a quantum computer by mapping the 
fermionic degrees of freedom to spin-1/2 particles [43 44|. Such a mapping is given by the famous Jordan- 
Wigner transformation. Here, the fermionic algebra can be expressed in terms of the su(2) algebra via c' k = 
— fef-Tj erf ) where = \ {uf, + The dynamical part of the fermionic many-body Hamiltonian 
often contains terms of the form hkj = c\.Cj + c^Ck, which become non-local after the transformation. 
Operators that are not adjacent in terms of the labeling often contain a chain of Pauli <j z operators in 
between them. A typical term of this kind that occurs after this transformation is = &k('8> :, iZk+i (7 i) a j- 
To simulate the time evolution of such a non-local term on a quantum computer, we need to be able to 
decompose this unitary into two qubit gates. Given the two unitaries Vki = cxp(ija^,af) and Ui = 
exp(ijaf) such a decomposition is indeed possible as depicted in the above circuit for the evolution of 
exp(-ieofcrfcrf). 
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